library(gstat)

# load up the required datasets
data(meuse)
data(meuse.grid)
data(meuse.riv)

#assign sensible shortcuts and manipulate them into sp objects
d <- meuse
class(d)
coordinates(d) <- ~ x+y
class(d)

grid <- meuse.grid
coordinates(grid) <- ~ x+y
gridded(grid) <- T
class(grid)

# this just illustrates that you should really look to other software
# for cartographic functionality
rv <- list("sp.polygons", SpatialPolygons(list(Polygons(list(Polygon(meuse.riv)),"meuse.riv"))), fill="lightblue")

#Let's look at the data
names(d)
summary(d$lead)
# lets look at lead

# plot the data
library(lattice)
trellis.par.set(sp.theme())
spplot(d, "lead", do.log=T, scales=list(draw=T))
bubble(d, "lead")

# normality?
densityplot(d$lead)
qqnorm(d$lead)
densityplot(log(d$lead))
qqnorm(log(d$lead))

# since we'll be dealing with this variable a bunch, make it stick
d$loglead <- log(d$lead)

# compute the sample variogram
v <- variogram(loglead ~ 1, data=d)
plot(v, plot.numbers=T)
# number of point needs to be > 100; override default lag width
v <- variogram(loglead ~ 1, data=d, width=130)
plot(v, plot.numbers=T)

# EXPLAIN 
# semivariance = measure of difference
# distance 
# nearby points are more similar

# Nugget, Sill, Range

# now we have to fit a model to the sample data
show.vgms()

# fit the model with parameters estimated "by hand"
vm1 <- vgm(psill=0.55, model="Sph", range=1000, nugget=0.05)
plot(v, model=vm1)

# A more objective method .. let gstat adjust the model to fit
vm2 <- fit.variogram(v, vm1)
vm2
plot(v, model=vm2)

# use that to predict the grid
g <- gstat(id="loglead", formula=loglead ~ 1, data=d, model=vm2)
okr <- predict(g, id="loglead", newdata=grid)

# graph it
spplot(okr, "loglead.pred", sp.layout=rv)

# write to raster dataset
library(rgdal)
writeGDAL(okr, "/home/perry/Desktop/meuse.tif", drivername="GTiff", type="Float32")

# write points to vector dataset
writeOGR(d, "/home/perry/Desktop/", "meuse", driver="ESRI Shapefile")

# open in Quantum GIS
q()
qgis /home/perry/Desktop/meuse.shp /home/perry/Desktop/meuse.tif

